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1  -  Numerical  Simulation  of  Viscous  Incompressible  Flows 


a  -  Continuation  methods 

Typically,  the  nonlinear  equations  resulting  from  discretizations  of  the  Navier-Stokes  equations  are  solved 
via  an  iterative  method  such  as  Newton’s  method  or  a  quasi-Newton  method.  For  interesting  values  of  the 
Reynolds  number,  these  iterative  methods  converge  only  when  a  sufficiently  good  initial  guess  is  available.  A 
popular  class  of  methods  for  generating  good  initial  guesses  are  continuation  methods  wherein  information 
obtained  at  lower  values  of  the  Reynolds  number  is  used  to  determine  the  initial  guess.  In  order  to  define 
a  continuation  method  one  must  select  a  step  size  in  the  Reynolds  number  and  then  choose  a  prediction 
algorithm  that,  using  previously  obtained  information,  determines  the  initial  guess. 

Our  research  has  shown  that  due  to  the  nature  of  solutions  of  the  Navier-Stokes  equations  it  essentially 
doesn’t  matter  which  prediction  algorithm  is  employed.  For  example,  one  can  take  the  same  step  size  in 
the  Reynolds  number  whether  one  simply  sets  the  inital  guess  to  be  the  solution  at  a  lower  value  of  the 
Reynolds  number  or  if  one  uses  a  more  complicated  tangent  plane  approximation  to  generate  the  initial 
guess.  Furthermore,  we  have  shown  that  the  allowable  step  in  the  Reynolds  number  is  roughly  proportional 
to  the  Reynolds  number  itself,  at  least  when  one  is  away  from  singular  points  of  the  solution. 

Our  results  are  based  on  some  new  estimates  obtained  for  the  path  derivatives,  i.e.,  derivatives  with 
respect  to  the  Reynolds  number,  of  solutions  of  the  Navier-Stokes  equations.  YVe  have  also  carried  out 
some  computational  experiments  that  support  our  analytical  results.  From  a  practical  point  of  view,  our 
results  indicate  that  it  is  most  efficient  to  use  the  simplest  type  of  predictor  algorithms  within  a  continuation 
method.  Furthermore,  a  rather  large  step  in  the  Reynolds  number  may  be  chosen. 

Details  of  the  above  results  may  be  found  in  the  first  item  in  the  bibliography  at  the  end  of  this  section. 

b  -  Parallel  algorithms  for  3-D  flows 

We  have  explored  the  feasibility  of  algorithms  for  three-dimensional  flows  that  can  take  advantage  of 
parallel  processing  architectures.  Our  algorithm  involves  the  choice  of  a  particular  form  of  the  governing 
equations,  a  finite  element  discretization,  and  an  advantageous  (from  a  parallel  processing  point  of  view) 
iterative  solution  strategy  for  solving  the  nonlinear  system  of  discrete  equations.  The  computational  problem 
is  reduced,  for  each  iteration,  to  solving  an  uncoupled  system  of  Poisson  equations.  Thus  we  have  some 
obvious  coarse-grain  parallelism  that  may  be  used  to  advantage  on  machines  (such  as  the  Cray-YMP)  that 
have  a  few  very  powerful  processors.  However,  since  one  must  solve  Poisson  equations,  one  can  also  take 
advantage  of  any  Poisson  solvers  having  fine-grain  parallel  structure  that  in  turn  can  make  good  use  of 
machines  with,  e.g.,  hypercube  architectures.  Among  candidate  algorithms  are  conjugate  gradient  methods 
and  multi-grid  methods. 

Our  own  efforts  so  far  have  been  directed  at  making  sure  that  the  combined  modeling,  discretization  and 
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solution  aspects  of  the  numerical  algorithm  result  in  practical  methods  from  the  points  of  view  of  accuracy, 
stability  and  efficiency.  This  was  a  necessary  exercise  since  many  facets  of  the  algorithms  we  consider  are 
“non-standard”  in  nature,  and  in  fact,  have  been  subject  to  considerable  debate  and  dispute  in  scientific 
circles.  Thus,  before  declaring  the  algorithms  to  be  viable  in  a  parallel  processing  environment,  we  had  to 
settle  some  of  these  outstanding  questions. 

The  algorithm  we  have  studied  is  a  velocity-vorticity  based  method.  For  two-dimensional  viscous 
incompressible  flow  problems,  the  streamfunction-vorticity  formulation  is  probably  the  method  of  choice 
for  use  in  conjunction  with  discretization  procedures.  For  three-dimensional  flows,  the  advantages  of  the 
streamfunction-vortivity  formulation  are  no  so  clear  cut.  Recently,  in  both  engineering  and  mathematical 
circles,  much  interst  has  been  focused  on  the  velocity-vorticity  formulation,  i.e. , 

div  u  =  0  , 


curl  u  =  L> 

and 

i/curlcurlw  +  curl(u>  x  u)  =  curlf , 

where  u  and  u  denote  to  unknown  velocity  and  vorticity  fields,  respectively,  and  u  and  f  are  the  given 
kinematic  viscosity  coefficient  and  the  body  force,  respectively.  The  central  problem  in  using  the  above 
formulation  is  determining  boundary  conditions  for  the  vorticity  at  boundaries  where  the  complete  velocity 
field  is  specified,  e.g.,  at  walls.  Note  that  the  above  equations  may  be  rewritten  in  the  form 

— Au  =  curl  w 

and 

—  A u  =  —curl (u  x  u)  +  curlf . 

We  have  carried  out  extensive  computational  experiments  using  the  velocity-vorticity  formulation.  Our 
studies  have  employed  finite  element  discretizations  of  three-dimensional  problems  and  have  focused  on 
question  of  accuracy  and  efficiency,  especially  in  view  of  different  choices  of  boundary  condition  treatment. 
Here  we  give  a  brief  account  of  what  we  have  learned. 

First,  whenever  the  velocity  is  given  at  a  boundary,  we  explicitly  the  information  given.  Some 

other  researchers  choose  to  only  use  the  normal  component  of  the  velocity,  leaving  the  tangential  boundary 
conditions  to  be  implicitly  satisfied.  Second,  for  the  vorticity,  the  best  thing  to  do  is  to  require  u  =  curlu 
at  boundaries  where  u  is  specified.  Note  that  u  •  n  =  curlu  ■  n  is  computable  on  a  boundary  where  u  is 
specified;  thus,  the  choice  we  make  is  to  set  the  tangential  components  of  the  vorticity  equal  to  the  tangential 
components  c'the  curl  of  the  velocity.  The  latter  requires  knowledge  of  interior  values  of  u.  We  have  found, 
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from  mathematical,  physical  and  computational  points  of  view,  that  the  above  combination  of  boundary 
conditions  are  the  prefered  ones. 

Our  next  result  is  concerned  with  the  accuracy  of  discrete  solutions.  We  have  found  that  in  order 
to  achieve  the  best  possible  accuracy  within  a  class  of  interior  discretizations,  one  must  approximate  the 
relation  u  =  curiu  to  better  accuracy  on  or  near  the  boundary.  This  may  be  accomplished  by  either  refining 
the  mesh  near  the  boundary,  or  by  using  higher  order  elements  at  the  boundary.  This  result,  i.e.,  that  one 
must  use  better  boundary  approximations,  has  not  been  previously  pointed  out  in  either  the  mathematical 
or  engineering  literature. 

Our  final  concern  with  the  velocity-vorticity  formulation  has  been  with  efficient  means  of  solving  the 
discrete  equations.  In  particular,  we  have  focused  on  iterative  schemes  that,  within  each  iteration,  uncouple 
the  velocity  and  vorticity  calculations.  Of  perhaps  greatest  interest  is  a  scheme  that  we  h»v«  found  to 
converge,  at  least  for  low  values  of  the  Reynolds  number,  and  which  requires  the  solution  of  six  uncoupled 
Poisson  equations  at  each  step  of  the  iteration. 

We  have  completed  a  paper  reporting  on  the  work  described  above;  the  paper  will  appear  as  an  invited 
chapter  in  the  fourth  volume  of  Computational  Methods  in  Viscous  Flows ;  it  is  the  second  paper  in  the 
references  listed  at  the  end  of  this  section. 

c  -  Magnetically  coupled  flows 

We  have  completed  the  analysis  of  finite  element  approximations  of  viscous  incompressible  flows  that 
are  affected  by  magnetic  and  electric  fields.  The  governing  equations  are  the  Navier-Stokes  equations  and 
the  Maxwell  equations,  fully  coupled  so  that  electromagnetic  fields  have  influence  on,  and  are  themselves 
influenced  by,  the  flow  field.  Along  the  way,  we  had  to  derive  new  existence  and  uniqueness  results  for  these 
coupled  equations.  One  feature  of  our  work  is  the  incorporation  of  a  variety  of  (inhomogeneous)  boundary 
conditions  into  both  the  mathematical  model  and  into  the  finite  element  methods  employed.  In  all  cases,  we 
have  been  able  to  prove  optimal  error  estimates  for  the  approximate  solution. 

At  the  present  time  we  aie  in  the  final  stages  of  running  computer  codes  that  have  been  developed. 
These  codes  not  only  illustrate  the  analytical  results  that  we  have  derived,  but  more  importantly,  have 
served  as  a  test  bed  for  devising  and  implementing  efficient  solution  algorithms  for  the  finite  dimensional 
nonlinear  system  of  equations  that  result  from  finite  element  discretizations  of  the  governing  system  of 
partial  differential  equations.  We  have  completed  the  preparation  of  a  paper  that  presents  the  analytical 
and  numerical  results  described  above;  this  manuscript  is  the  third  one  in  the  list  of  references  at  the  end  of 
this  section. 
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2  -  Numerical  Approximation  of  Control  Problems 
a  -  Control  of  viscous  incompressible  flows 

We  have  completed  the  analysis  of  various  theoretical  and  computational  issues  associated  with  the 
control  of  viscous  incompressible  flows,  specifically  those  described  by  the  Navier-Stokes  equations.  All 
types  of  controls  have  been  considered.  From  the  point  of  view  of  analysis,  the  easiest  cases  to  treat  were 
distributed  controls,  e.g.,  body  force  controls,  and  Neumann  controls,  e.g.,  stress  controls.  However,  the 
more  interesting  case  is  that  of  Dirichlet  controls,  e.g.,  control  by  adjusting  the  velocity  at  the  boundary. 
The  latter  has  application  in,  e.g.,  drag  control  by  blowing  or  sucking  fluid  through  the  boundary.  In  all 
cases,  one  of  our  aims  was  to  minimize  drag  by  appropriately  choosing  the  control. 

Typical  of  the  generality  of  the  problems  we  have  treated  is  the  following  situation.  We  have  that  the 
state  variables  and  controls  must  satisfy  the  Navier-Stokes  system 

-t/div  ((grad  u)  +  (grad  u)T )  +  u  •  grad  u  +  grad  p  =  f  +  gd  in  Q , 

div  u  =  0  in  Q  , 
u  =  ho  4-  go  on  T! 

and 

~pn  +  i/(grad  u  +  grad  uT)  •  n  +  i/au  =  hN  +  gN  on  T2  . 

Here  is  a  bounded  two  or  three-dimensional  domain  and  Tj  and  r2  are  parts  of  its  boundary.  The  state 
variables  are  the  velocity  u  and  the  pressure  p.  The  functions  f.a.hp  and  hjv  are  given.  The  controls  are 
denoted  by  gd  (distributed  control),  go  (Dirichlet  control)  and  g^v  (Neumann  control).  The  controls  need 
not  act  on  all  of  their  respective  domains;  e.g.,  gd  may  act  on  only  a  subset  of  fi  and  go  may  act  on  only  a 
subset  of  IV  Our  optimization  problem  can  involve  choosing  controls  such  that  either 

<7(u)  =  \  J  lu  “  u0|2dfi 

or 

£(u)  =  ^  J  Kgrad  u)  +  (grad  u)T  |2  <fD 

is  minimized.  In  the  first  case  we  want  the  velocity  field  to  be  as  “close”  as  possible  to  some  prescribed  field 
uo;  in  the  second  case  we  are  minimizing  the  viscous  drag. 

For  the  case  of  Dirichlet  controls,  we  have  treated  controls  that  are  unknown  functions  defined  un 
portions  of  the  boundary  of  the  flow  region.  We  have  also  treated  problems  where  the  fluid  is  allowed  to 
enter  or  leave  the  flow  region  only  through  a  finite  number  of  specified  holes  on  the  boundary,  and  with 
a  specified  velocity  distribution  through  the  holes.  The  only  control  is  the  amount  of  fluid  that  can  enter 
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or  exit  through  each  hole.  This  last  problem  effects  considerable  simplifications  since  the  set  of  admissible 
controls  is  finite  dimensional. 

We  have  fully  analyzed  questions  of  existence  of  optimal  controls  and  of  Lagrange  multipliers  for  all 
cases.  We  have  also  developed  and  analyzed  finite  element  algorithms;  all  the  agorithms  we  have  developed 
yield  optimally  accurate  approximations.  At  the  present  time  we  are  implementing  some  of  these  algorithms 
into  a  computer  code  with  our  main  goal  being  developing  efficient  solution  algorithms  for  the  formidable 
discrete  problem  that  results  from  the  discretization  of  the  control  problem. 

This  work  is  been  reported  in  papers  4-7  listed  at  the  end  of  this  section. 

b  -  Boundary  controllability  problems  for  the  wave  equation 

We  have  been  trying  to  develop  algorithms  for  approximately  solving  controllability  problems  for  the 
wave  equation.  Our  eventual  goal  is  to  extend  these  algorithms  to  controllability  problems  for  elastic  and 
viscoelastic  structures,  including  large  space  structures. 

For  problems  in  one  space  dimensnion,  one  can  easily  devise  numerous  algorithms  that  can  be  used  to 
approximately  solve  controllability  problems  fo  the  wave  equation.  In  higher  dimensions,  the  task  is  not  so 
easy.  We  believe  we  have  developed  and  effective  means  for  treating  such  controllability  problems.  Consider 
the  model  problem 

u„  -  A u  =  0  in  (0,T)  x  fl 
u(0,x)  =  u0(x)  for  x  € 
u<(0,  x)  =  «i(x)  for  x  6 
u(T,  x)  =  uz(x)  for  x  6  fl 
u,(T,  x)  =  U3(x)  for  x  e  ft  ■ 

Here  Uj,  j  =  0,...,3,  are  prescribed  functions.  It  is  well  known  that  for  sufficiently  large  T  there  exists 
boundary  controls,  either  of  Neumann  or  Dirichlet  type,  such  that  the  above  problem  has  a  solution.  In  fact, 
in  general,  the  solution  is  not  unique. 

Our  approach  is  to  solve  the  above  underdelerrnined  problem,  or  rather  discrete  versions  of  the  above 
problem,  by  a  least  squares  techniques.  This  results  in  sparse  positive  definite  discrete  problems  to  be  solved 
for  which  efficient  algorithms  are  known.  At  the  present  time  we  are  developing  codes  that  implement  these 
methods  with  a  view  towards  choosing  the  best  possible  among  variants  of  the  basic  least  squares  procedure. 

A  preliminary  report  describing  this  work  has  been  completed;  it  is  paper  number  8  in  the  list  of 
references  at  the  end  of  this  section. 

3  -  Centroidal  Voronoi  Grids 

Voronoi  grids  and  their  duals,  the  Delaunay  triangulations,  are  now  in  widespread  use  in  the  numerical 


computation  of  a  variety  of  problems,  ranging  from  partial  differential  equations  to  image  processing  to  data 
compression.  These  tesselations  are  simple  to  define.  Given  a  set  of  points  z j  &  Cl  C  fftn,  the  Vorornoi 
region  for  a  point  z*  is  the  set  of  points  in  Cl  that  are  closer  to  z*  than  to  any  of  the  other  given  points. 
Voronoi  regions  turn  out  to  be  polygons  and  their  usefulness  in  the  numerical  solution  of  partial  differential 
equations  stems  from  the  fact  that  the  dual  grid  is  in  general  a  triangulation  of  the  given  points  z;  with  the 
property  that  among  all  possible  triangulations  of  these  points,  it  is  the  one  with  the  maximum  minimum 
angle.  This,  of  course,  is  a  desirable  property.  Furthermore,  there  has  been  much  research  and  success 
recently  into  developing  efficient  algorithms  for  computing  Voronoi  grids. 

In  general,  the  centroids  y ,  of  the  Voronoi  polygons  need  not  be  anywhere  near  the  points  Zj  from 
which  the  Voronoi  polygons  are  determined.  However,  we  want  to  consider  the  case  where  y j  =  z j.  More 
generally,  given  a  region  Cl  €  Rn,  a  positive  integer  AT,  and  a  density  function  p(x),  we  want  to  find  TV 
points  Zj ,  j  =  1, . . . ,  TV  and  TV  non-overlapping  regions  Clj,  j  =  1  covering  Cl  such  that  simultaneously 

each  Qj  is  a  Voronoi  region  for  zj,  and  each  z,  is  a  centroid,  with  respect  to  p,  of  Clj.  Thus,  the  set 
{zj,Clj},j  =  1, . . . , TV,  should  satisfy 

|x  —  Zj  |  <  |x  —  z*  |  Vx  €  Clj  and  all  k  j 


and 


-  in,  dQ 
Z'~  In 


for  j  =  1, . . . ,  TV  . 


We  have  been  able  to  show  that  this  problem  has  a  (non-unique)  solution.  Furthermore,  we  have  also 
considered  the  following  algorithm,  known  as  Lloyd’s  algorithm  in  the  literature,  and  shown  that  it  produces 
a  sequence  of  iterates  that  converge  to  a  solution.  We  begin  with  an  arbitratry  set  of  points  zj  £  that 
are  in  general  position.  At  the  m-th  step,  we  start  with  the  points  zj*-1  and  then  compute  the  Voronoi 
tesselation  of  Q  for  these  points,  which  we  denote  by  Cl™.  This  step  can  be  acomplished  by  well  known 
existing  algorithms.  We  then  define  z™  to  be  the  center  of  mass,  with  respect  to  the  given  density  function 
p(x),  of  the  region  Cl™.  We  have  shown  that  this  algorithm  converges  in  general  and  have  applied  it  to  some 
practical  data  compression  and  image  processing  problems.  We  have  also  studied  a  variant  of  the  above 
algorithm  wherein  the  explicit  construction  of  the  Voronoi  regions  is  avoided  in  favor  of  a  stochastic  process 
that  yields  the  same  result. 

We  have  also  been  interested  in  numerical  analysis  applications  of  centroidal  Voronoi  grids.  For  example, 
here  is  one  such  application.  Consider  a  function  /(x)  and  its  weighted  integral  over  a  given  region  Cl.  Now 
consider  piecewise  constant  type  quadrature  rules  of  the  type 


where  Zj  €  Q;  and  Vj,j  =  1, . . . , N,  are  the  volumes  of  the  regions  which  tesselate  Q.  Then,  for  the  class 
of  Lipschitz  continuous  functions,  among  all  possible  choices  of  quadrature  regions  Qj  and  qudarature  points 
Zj ,  the  one  that  yields  the  most  accurate  quadrature  rule  is  one  for  which  the  z's  and  Q's  form  a  centroidal 
Voronoi  grid.  Another  application  of  centroidal  Voronoi  grids  we  have  examined  is  to  the  definition  of  finite 
difference  schemes  on  general  grids.  It  turns  out  that  easily  defined  finite  difference  schemes  that  are  only 
first  order  accurate  on  arbitrary  grids,  are  second  order  accurate  on  centroidal  Voronoi  grids;  since  the  latter 
are  perfectly  useful  for  tesselating  arbitrary  domains  ft,  we  see  that  centroidal  Voronoi  grids  are  useful  in 
generating  finite  difference  schemes  on  arbitrary  domains. 

We  are  currently  writing  a  paper  presenting  these  results;  it  is  paper  9  in  the  reference  list  at  the  end 
of  this  section. 

4  -  A  Book  on  Finite  Element  Methods  for  the  Navier-Stokes  Equations 

We  have  completed  the  writing  of  a  book  on  finite  element  methods  for  the  Navier-Stokes  equations. 
This  book  is  an  outgrowth  of  much  of  the  research  we  have  carried  out  under  AFOSR  sponsorship.  Thus,  it 
draws  form  our  own  considerable  expereince,  as  well  as  that  or  others.  The  goal  of  the  book  is  to  acquaint  the 
interested  reader,  be  it  an  engineer  or  mathematician,  with  the  mathematical  results  and  the  best  algorithms 
available  for  solving  the  Navier-Stokes  equations  by  finite  element  methods.  The  book  appeared  at  the  end 
of  the  summer  of  1989;  the  exact  citation  is  found  below  at  the  end  of  the  list  of  references. 
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